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Abstract 

The theory underlying a proposed random number generator for numeri- 
cal simulations in elementary particle physics and statistical mechanics is dis- 
cussed. The generator is based on an algorithm introduced by Marsaglia and 
Zaman, with an important added feature leading to demonstrably good sta- 
tistical properties. It can be implemented exactly on any computer complying 
with the IEEE-754 standard for single precision floating point arithmetic. 
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1. Introduction 



Numerical simulations in elementary particle physics and statistical me- 
chanics are increasingly performed on massively parallel computers. These 
machines offer unmatched computing power, thus making it possible to simu- 
late larger systems and to achieve greater statistical precision. It is well-known 
that the random number generators employed in these computations can be a 
source of systematic error. In fact many of the popular generators used to date 
failed to give correct results in some recent simulations of the 2-dimensional 
Ising model [1,2]. While the Ising model is a rather special case, with un- 
usual regularity, the lesson clearly is that random number generators should 
be chosen with care, especially when one aims for high-precision results. 

The generator discussed in this paper derives from an algorithm originally 
proposed by Marsaglia and Zaman [3]. It has a very long period and excellent 
statistical properties on short and long time scales. The quality of the gen- 
erator is established using some mathematical results on chaotic dynamical 
systems, the spectral test and a number of empirical tests. 

The algorithm has been implemented on the APE-100, a parallel computer 
now intensively used in elementary particle physics (for a short description and 
guide to the literature see ref.[4]). One may also easily write a FORTRAN 
code for the generator, which will run correctly on any machine complying 
with the IEEE-754 standard for single precision floating point arithmetic. 

The definition and basic properties of the Marsaglia- Zaman algorithm 
are reviewed in sect. 2. For appropriately chosen parameters the period of 
the generator can be proved to be very large [3]. Its statistical properties 
are however not as good as initially assumed. In particular, the generator 
fails in the classical gap test [5] and an unfavourable lattice structure in the 
distribution of random numbers in high dimensions has been discovered [6,7]. 

The important new observation made in this paper is that the Marsaglia- 
Zaman algorithm is closely related to a dynamical system, which is known to 
be chaotic in a strong sense (it is a so-called li'-system [11]). One then infers 
that the correlations detected in the gap test, for example, are short ranged in 
time. A sequence of random numbers with much better statistical properties 
is therefore obtained by picking out elements of the original sequence at time 
intervals greater than the correlation time. All this is explained in sects. 3 and 
4, and the quality of the so improved generator is evaluated in sect. 5. Imple- 
mentation details and timing benchmarks for various machines are included 
for completeness. 
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2. Marsaglia-Zaman generator 



The random number generator defined below is based on a so-called 
subtract-with-borrow algorithm [3]. For the particular choice of parameters 
specified in subsect. 2.4 the generator is known by the name of RCARRY [8]. 

2.1 Definition 

Let b be an arbitrary integer greater than 1, referred to as the base, and define 
X to be the set of integers x satisfying < x < b. The algorithm generates 
a random sequence xo, x\, X2, ■ ■ ■ of elements of X recursively, together with 
a sequence Co, ci, C2, . . . of "carry bits". The latter take values or 1 and are 
used internally, i.e. the interesting output of the algorithm are the numbers x n , 
or rather x n /b, if one requires random numbers uniformly distributed between 
and 1. 

The recursion involves two fixed lags, r and s, which are assumed to 
satisfy r > s > 1. For n > r one first computes the difference 

^■n — %n — s %n — r ^n — 1 ? (2T) 

and then determines x n and c n through 

x n = A n , c n = if A n > 0, 

(2.2) 

x n = A n + b, c n = 1 if A n < 0. 

It is trivial to verify that x n is contained in X if x n _ s and x n _ r are and if 
c n _i is or 1. The name "carry bit" for c n is now quite intuitive, since c n 
simply indicates whether a shift by the base b was necessary when computing 
x n . 

To start the recursion, the first r values xo, x\, . . . , x r _\ together with an 
initial carry bit c r _i must be provided. The configurations 

x = x\ = . . . = x r _i = 0, c r _i = 0, (2-3) 

x = x\ = . . . = x r _i = 6 — 1, c r _i = 1, (2-4) 

should be avoided, because the algorithm yields uninteresting sequences of 
numbers in these cases. All other choices of initial values are admitted in the 
following and we shall then say that the generator has been properly initialized. 
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2.2 Period of the generator 

For some values of the base b and the lags r,s, the period of the sequence 
generated through eqs.(2.1),(2.2) can be determined rigorously. Define 



to 



= ¥ - 



V + 1 



(2.5) 



and let q be the smallest positive integer such that 



b q = 



1 mod to. 



(2.6) 



The existence of q is guaranteed since to and b are relatively prime. An 
important mathematical result of Marsaglia and Zaman now is [3] 

Theorem 2.1. If to is a prime number, the period of the generator dehned 
through eqs.(2.1),(2.2) is equal to q. More precisely, if the generator has been 
properly initialized, the following is true. 

1. For all n > r we have x n+q = x n . 

2. Any number p, such that x n+p = x n for more than r successive values of 
n, is an integer multiple of q. 

It should be emphasized that the period is independent of the chosen initial 
values xo,x\, . . .,x r _\. Note that this particular string of numbers may not 
occur anywhere else in the sequence, i.e. in general the algorithm gets into a 
loop only after the first r updates have been made. 

Another comment is that the period of the generator must be expected 
to depend on the initial values, if to is not prime. Such generators are not safe 
and should be avoided unless all periods can be shown to be large. 

2.3 Associated linear congruential generator 

The algorithm of Marsaglia and Zaman is closely related to the standard linear 
congruential generator with multiplier 



and modulus to [6]. Such generators have been studied vigorously in the past 
and we shall later rely on some of this theory when we discuss the statistical 
properties of the random number sequence produced by the Marsaglia- Zaman 
algorithm. 





3 



The linear congruential generator alluded to above operates on the set of 
all integers y in the range < y < m. Starting from an initial value yo, a 
sequence of random numbers yo, y\, y2, ■ ■ ■ is obtained recursively through 

y n = ay n -\ mod to. (2.8) 

The multiplier a satisfies 

ab = 1 mod to (2-9) 
and the recursion is thus equivalent to 

y n = by n+1 mod to. (2.10) 

It is not difficult to show that the period of the sequence is equal to q if to is 
prime. 

The relation between this generator and the Marsaglia-Zaman generator 
is summarized by [6] 

Theorem 2.2. Let (x n ) n >o be a sequence of random numbers generated 
through the Marsaglia-Zaman algorithm, with carry bits (c n ) n >o and proper 
initial values. Then, for all n > r, the integers 

T—l S—l 

Vn = x n _ r+k b k - ^2 x n - s +kb k + c n _i (2.11) 

k=0 k=0 

are in the range < y n < m. Moreover the relation 

by n+1 - y n = mx n (2.12) 

holds and the sequence (y n ) n > r is thus generated through the recursion (2.8). 

The theorem shows at once that the Marsaglia-Zaman algorithm is essentially 
a clever way to implement certain linear congruential generators with huge 
moduli. Manipulations of large integers are avoided by breaking them up into 
a vector of smaller numbers which are then processed one by one. 
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2.4 Choice of parameters 

Most computers used for large scale numerical simulations have been designed 
to yield maximum performance for floating point operations. The parameters 
b, r and s should thus be chosen so as to be able to implement the generator 
using floating point arithmetic. 

Single precision real numbers on computers complying with the IEEE- 
754 standard are represented by a string of 32 bits, with 23 bits reserved for 
the mantissa and the rest for the sign and exponent of the number. Signed 
integers of absolute magnitude up to 2 24 can thus be dealt with exactly on 
such machines using floating point arithmetic. So if we choose 

6 = 2 24 , (2.13) 

all elements of X (and b itself) will be computer representable numbers. As 
for the lags r and s, we take 

r = 24, s = 10, (2.14) 

a choice proposed by Marsaglia and Zaman [3] and recommended by James 
[8]. The difference A n in the recursion (2.2) then is 

A n = X n _ W - X n _2A ~ Cn-l, (2.15) 

and 24 integers xo,x\, . . .,#23 m the range < Xk < 2 24 plus a carry bit C23 
are required to initialize the generator t. Note that no rounding occurs in the 
computation of A n , since the final and intermediate results are representable 
numbers, i.e. the algorithm is implemented exactly. 

The modulus m and multiplier a for this choice of parameters are given 

by 

m = 2 576 - 2 240 + 1, (2.16) 
a = 2 576 - 2 552 - 2 240 + 2 216 + 1. (2.17) 

Using elementary number theory and the complete decomposition of m — 1 
into prime factors, it is possible to prove that m is a prime number [3]. The 

I The FORTRAN code for this algorithm printed in ref.[8] contains an error. A correct 
program is obtained by interchanging the indices 124 and J24 in the line UNI=SEEDS(I24)- 
SEEDS(J24)-CARRY [9]. 
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period of the generator is thus determined by theorem 2.1. Some further work 
then yields 

q = (m- l)/48 ~ 5.2 X 10 m , (2.18) 

which is a very long period indeed. There is no chance that, on any earthly 
computer, one will ever come close to exhausting this sequence of random 
numbers. 

In the following the parameters of the generator are assumed to be as 
specified above. The reader should however meet no difficulty in carrying over 
the discussion to any other case of interest. 



3. Origin of statistical correlations 

The Marsaglia-Zaman generator is now known to fail in several empirical 
tests of randomness, a particularly simple case being the gap test ([5]; for a 
lucid description of the test see ref.[10], p.60f). As explained below there are in 
fact some rather obvious correlations between successive vectors of r random 
numbers. They are seen most clearly when the generator is described in the 
language of dynamical systems. 

3.1 Geometrical preliminaries 

The unit hyper-cube in r dimensions is the set of all vectors 

v = (v ,v 1 , . . .,w r _i) (3.1) 

with real components between and 1. If opposite faces of the hyper-cube are 
identified one obtains an r dimensional torus T' r . The points on this manifold 
are also represented by vectors v, as above, with the understanding that v and 
w describe the same point if v k = w~k mod 1 for all h. 

T r contains a discrete subset, T r , which consists of all vectors v with 
components of the form 

Vk = n k /b, n k = 0, 1,2, . . .,b - 1. (3.2) 

T r is an r dimensional hyper-cubic lattice with spacing 1/6, which may be 
regarded as a discrete approximation of the torus. 
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The distance between any two points v and w on T' r is denned through 
d(v,w) = maxd k , d k = mm{\v k - w k \, 1 - \v k - w k \\. (3.3) 

k 

It is straightforward to check that d has all the properties required for a 
decent distance function on T r . In particular, it is invariant under translations 
modulo 1. 

3.2 The Marsaglia-Zaman generator as a dynamical system 

Let us now consider a sequence of random numbers xo, x\, X2, ■ ■ ■ generated 
through the Marsaglia-Zaman algorithm, with carry bits (c n ) n >o and proper 
initial values. The vectors 

v(t) = (x n ,x n+1 , . . .,a; n+r _i)/6, n = rt, (3.4) 

define a point on the (discrete) torus T r which moves as the "time" / progresses 
from in steps of 1. If we also introduce a time dependent carry bit, 

c(t) = C rt+r -l, (3.5) 

it is clear that the evolution of v(t) and c(t) is determined by the recursion 
(2.1),(2.2). 

We are thus led to interpret the Marsaglia-Zaman generator as a discrete 
dynamical system, consisting of a set S of states and a mapping (f> : S S . A 
state is defined by a point on the discrete torus and a carry bit. (f> maps any 
such state onto the next one, viz. 

(v(t+l),c(t+l))=<t>{v(t),c(t)). (3.6) 

Note that (f> does not refer to any of the previous states. One only needs to 
know the current state to be able to compute the next one. 

3.3 Continuity and statistical correlations 

For a good generator one requires that successive vectors of random numbers 
be statistically independent. That is, if (v,c) runs through all possible states, 
the joint distribution of (v , c) and (f>(v, c) should be uniform on S X S . 

Of course this cannot be true since (f> operates on a finite set of states. The 
distribution is at best approximately uniform. Since one can only generate a 
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relatively small number of states in practice, one is anyway unable to test the 
distribution very precisely. One should however be worried by correlations 
that are strong enough to give a measurable effect in any simple statistical 
test. 

We now show that such correlations exist. Let us first ignore the carry 
bits. The recursion (2.1), (2. 2) then reads 

x n = x n _ s - x n _ r mod b (3.7) 

and (f> becomes a linear transformation of the torus. An important consequence 
of this fact is that nearby points are mapped onto nearby ones. So if one 
chooses a set of random points v in some small volume, their successors (f>(v) 
are contained in some other small volume. In particular, they are not scattered 
over the whole torus, as one would expect if (f>(v) were statistically independent 
of v. 

The carry bits only affect the least significant digits of the random num- 
bers and so cannot destroy the basic continuity of (f>. More precisely, if we 
define 

(v,c) = cf>(v,c), (3.8) 

it is possible to show that 

d(v, w) < 4d(v, w) + 3/6. (3.9) 

The distance between two points on T' r thus increases by at most a factor 4 
plus 3 lattice spacings. In particular, small regions are mapped onto small 
regions and so we again conclude that successive vectors of random numbers 
are strongly correlated. 

It should be emphasized that the effects caused by these correlations are 
readily seen in empirical tests. In particular, the failure of the Marsaglia- 
Zaman generator in the gap test can be explained in this way. Note, inciden- 
tally, that similar correlations are present in all lagged Fibonacci generators 
using addition or subtraction as the binary operation. 
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4. Deterministic chaos 



A characteristic feature of chaotic dynamical systems is that trajectories 
starting at nearby states diverge exponentially with time. Even if the evolution 
is locally continuous, such a system appears to behave randomly on larger time 
scales. One could also say that any state specified to some finite precision has 
an exponentially deteriorating memory of its history. We now show that the 
dynamical system underlying the Marsaglia-Zaman generator is chaotic in this 
sense. 

4-1 Numerical experiment 

It is helpful to start with a simple experiment illustrating the chaotic nature of 
the mapping (f>. The experiment consists in choosing a random sample of 1000 
pairs of trajectories (v(/),c(/)) and (v '(/), c'(/)) , with initial values separated 
by 1 lattice spacing, viz. 

d(v(0),v'(0))=l/b. (4.1) 

One then computes the average distance 

S(t) = {d(v(t),v'(t))) (4.2) 

as a function of the evolution time /. 

Fig. 1 shows that the trajectories are rapidly diverging. In the range 
4 < t < 16 the data are well described by 

6(t) = Ae\ A = 5 X 10" 8 , (4.3) 

i.e. the separation is growing exponentially with a rate close to 1. Around 
/ = 17, S(t) levels off and assumes a value equal to 12/25 within statistical 
errors. This is the average distance between two randomly chosen points on 
the torus, thus indicating that v(t) and v'(t) are no longer correlated. 
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Fig. 1. Average distance 6(t) between neighbouring trajectories as a 
function of the evolution time t. 

4-2 Continuum limit 

For the further study of deterministic chaos it is now useful to pass to the 
continuum limit 1/6 — > 0, where the space of states S becomes equal to the 
full torus T' r and the carry bit is neglected. This is an accurate approximation 
to the discrete system on short time scales and if all distances of interest are 
much greater than the lattice spacing. In particular, the evolution of diverging 
trajectories can be expected to be correctly described when they are sufficiently 
far apart. 

In the continuum limit the mapping (f> reduces to 

<f)(v) = L r v mod 1, (4.4) 
where L is the linear transformation 

Lv = Oi,-u 2 , • • . ,w r _i,?V_ s - v ). (4.5) 
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L can be considered an r X r matrix with entries 0, 1 and —1. It is then trivial 
to verify that det L = 1 and (f> is hence invertible and volume preserving. 

According to the established mathematical terminology, the continuum 
system (T r ,/2,(j)) (where fj, denotes the standard measure on T r ) is a classical 
dynamical system. The occurence of chaos in such systems has been studied 
extensively and many deep results have been obtained. In the rest of this 
section the system (T' r ,/2,(j)) will be discussed from the point of view of the 
mathematical theory. Although no previous knowledge on dynamical systems 
is required, the reader may now find it useful to consult one or the other book 
on the subject such as refs. [11-13], for example. 

4-3 Liapunov exponent 

In the continuum system the exponential rate of divergence of neighbouring 
trajectories can be computed analytically as follows. 

Suppose v(t) and v'(t) are two trajectories such that their distance is very 
much smaller than 1 at / = 0. Let us define the difference vector 

u(t) = v'(t) - v(t) mod 1, ~\<u k {t)<\. (4.6) 

It is clear that the norm of this vector, 

||ti(i)|| = max|ti fc (i)|, (4.7) 

k 

is equal to the distance between the trajectories at time /. Furthermore, from 
eq.(4.4) one infers that 

u{t + 1) = L r u{t) (4.8) 

if < a condition which is satisfied as long as the trajectories are 

sufficiently close. 

The dominant exponential growth of the deviation vector u(t) is hence 
determined by the largest eigenvalues of L. The characteristic equation of L, 

A r -A r " s + 1 = 0, (4.9) 

can easily be solved numerically and one finds that all eigenvalues A are com- 
plex and non-degenerate. There are 4 eigenvalues with maximal absolute value 
given by 

|A| max = 1.04299... (4.10) 
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Now if the initial deviation vector m(0) has a non-zero component in the di- 
rection of the corresponding eigenvectors (which is the generic case), one con- 
cludes that 

||ti(i)|| oc e vt (4.11) 

at large times /, where 

z/ = rln|A| max = 1.01027... (4.12) 

Of course eq.(4.11) only holds as long as the evolution equation (4.8) applies. 
By considering smaller and smaller initial deviations, this condition will be 
fulfilled for any desired length of time. Eq.(4.11) then becomes asymptotically 
exact. 

The exponent v is referred to as the Liapunov exponent of the system. 
As already noted in subsect. 4.2, the evolution of diverging trajectories in 
the discrete system is expected to be accurately described by the continuum 
system. A comparison of the result of the experiment, eq.(4.3), with the value 
of the Liapunov exponent v confirms this. We have thus shown that the 
chaotic behaviour of the Marsaglia-Zaman generator can be traced back to 
the instability of the underlying lagged Fibonacci generator. 

4-4 Kolmogorov entropy and mixing 

The continuum system (T r , fj,, (f>) can be proved to belong to a class of strongly 
unstable systems. While the relevance of this remark for the discrete system is 
not completely obvious, it does provide some further insight into how repeated 
application of a smooth mapping can lead to randomness. 

The mapping (f> is in many respects similar to the famous cat map of 
Arnold. In particular, under the action of (f> the torus is stretched in r/2 
directions and shrunk in r/2 complementary directions. After many iterations 
any region in T' r (a cat's body, for example) is first made very long and thin 
and then wrapped on the torus. As a result the region is scattered over the 
whole manifold. 

These heuristic remarks can be made much more precise and it is then 
possible to show, using the theorems discussed in ref.fll], that (T r ,/2,(j)) is a 
so-called li'-system. This means that it has a positive Kolmogorov entropy 
and that consequently it is mixing and ergodic. 

The property of mixing is particularly intuitive. It states that 

lim n (AC) <j>\B)) = fj,(A)fji(B) (4.13) 

t^oo 
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for all measurable sets A,B. In other words, if the set B is evolved for a 
long time, it will be uniformly distributed over the torus and thus occupies a 
fraction fJ,(B) of every other set A (recall that (f> is volume preserving). 

The Kolmogorov entropy is a substantially more difficult notion. Basically 
it is the rate at which the knowledge about the system is lost as it evolves from 
an only imprecisely specified initial state. A positive entropy thus implies that 
one loses information exponentially fast. 



5. Improved generator 

The important qualitative implication of the chaotic nature of (f> is that 
the correlations discovered in sect. 3 are short ranged in time. A sequence 
of random numbers with significantly better statistical properties is therefore 
obtained by keeping only a fraction of the full sequence of numbers produced 
by the Marsaglia-Zaman algorithm. The precise rule is given below and several 
statistical tests are performed to confirm the expected improvement. 

5.1 Definition 

We again start from a sequence of random numbers xo, x\, X2, ■ ■ ■ generated 
through the Marsaglia-Zaman algorithm, with carry bits (c n ) n >o and proper 
initial values. Instead of using all numbers x n , we now read r successive 
elements of the sequence, discard the next p — r numbers, read r numbers, and 
so on. The integer p > r is a fixed parameter which allows us to monitor the 
fraction of random numbers "thrown away". In particular, the old generator 
corresponds to p = r, where no numbers are discarded. 

The numbers selected in this manner define a history of states [v(t), c(/)) 
through 

v(t) = (x n , x n+1 ,. . . , s n+r _i)/6, n = pt, 

(5.1) 

c(t) = C n+r _i. 

As before the time evolution is generated by a well-defined mapping <f), p : S ^ S 
such that 

(v(t+l),c(t+l)) =<f> p (v(t),c(t)). (5.2) 
In the continuum limit (f> p reduces to the linear transformation 

<f> p (v) = L p v mod 1, (5.3) 
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where L is given by eq.(4.5). 

The discussion in sect. 4 now suggests that deterministic chaos leads to a 
complete decorrelation of successive states for values of p greater than about 
16r = 384. For such p the corresponding sequence of random numbers is 
expected to possess excellent statistical properties. In practice one may be 
satisfied with a smaller value of p, as a full decorrelation, down to the level of 
the least significant bits, may in many cases be unnecessary. The statistical 
tests reported in the following subsections help to clarify the situation and a 
more definite recommendation as to which value of p to choose will be issued 
after that. 

5.2 Spectral test 

For any state (v,c) an integer y in the range < y < to may be defined 
through 

r— 1 s— 1 

y = J>*& fc+1 -Y^vr-s+kb^ 1 + c, (5.4) 

fc=0 fc=0 

where vo, v\, . . . , f r -i are the components of v (cf. theorem 2.2). y should be 
regarded as an observable constructed from the given state. In particular, a 
trajectory (v(/),c(/)) of states, generated by the mapping <f> v , is associated 
with a sequence of values y(t). Theorem 2.2 tells us that 

y(t + 1) = a p y(t) mod to, (5-5) 

i.e. (f> p is related to a linear congruential generator with modulus to and mul- 
tiplier a' p mod to. 

The multi- dimensional distributions of y(t) can be studied by applying the 
powerful spectral test for linear congruential generators. The test effectively 
probes the statistical independence of successive states c(/)), since any 

correlation between the values of y(t) can be regarded as a correlation among 
the corresponding states. For a detailed description of the spectral test the 
reader is referred to Knuth's book [10]. Here we merely introduce the necessary 
notations and discuss the results of the test. 

An infamous property of linear congruential generators is that vectors of 
D successive random numbers fall on parallel hyper-planes with often appre- 
ciable spacing. The spectral test consists in calculating the maximal spacing 
ho, or rather the "accuracy" = l/hp, for low dimensionalities D. The 
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Table 1. Merits /ic of some generators with modulus m and multiplier aP mod m 



p 


M2 


M3 


M4 


M5 


Me 


M7 


Ms 


24 


4e 29 


e 85 


e 56 


2e 27 


e 86 


3e 72 


g e 58 


48 


0.20 


0.07 


0.03 


9e 23 


5.08 


2e 33 


2e 31 


96 


2.67 


1.04 


1.64 


0.04 


1.60 


0.14 


0.10 


192 


1.82 


0.67 


0.70 


1.53 


2.69 


4.78 


1.54 


384 


0.56 


0.82 


2.30 


1.56 


0.84 


4.60 


0.29 


768 


1.63 


2.59 


3.08 


0.59 


0.96 


1.29 


1.12 


223 


1.80 


0.87 


2.39 


3.79 


2.29 


0.78 


2.29 


389 


2.27 


3.46 


3.92 


2.49 


2.98 


4.23 


0.46 




m and a are 


given by eqs.(2.16),(2.17)] 









outcome of the spectral test may be rated through the figures of merit 

M-D = r n n , -iv 5 - 6 
ml {jL> + 1) 

Good generators achieve values of fj,£> greater than 1 for say D = 2, . . . , 6. On 
the other hand, if the merit is significantly smaller than 0.1 for some of these 
dimensions, one has picked a particularly bad multiplier. 

The results of the spectral test are listed in table 1. The first line cor- 
responds to the original generator where no random numbers are discarded. 
As already noted in refs.[6,7], there are strong correlations between successive 
values of the observable y in this case, for any dimensionality D. Evidently 
this generator is a poor source of random numbers. 

In general the merits are quite acceptable for p greater than about 200. 
The merits for two favoured values around 200 and 400 are listed in the last 
two lines of table 2. All this is very much in line with what one expects from 
deterministic chaos. It should however be emphasized that the spectral test is 
a full period test, while the decorrelation through diverging trajectories takes 
place on short time scales. 
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5.3 Further statistical tests 

a. Serial correlation test. This test is applied to the associated linear con- 
gruential generator. It is a full-period theoretical test, where one computes 
the correlation coefficient between successive values of y exactly (see ref.[10] 
for further explanations). For values of p greater than about 100 it is passed 
easily. 

b. Gap test. In ref.[5] the original generator (p = 24) has been subjected to a 
large number of empirical tests. All tests were passed with the exception of 
the gap test. This test has now been repeated for various values of p, with the 
same test parameters, and no significant statistical correlations were detected 
for p > 48. 

c. Ising model. Simulations of the 2-dimensional Ising model, using cluster 
algorithms, have proved to be a particularly sensitive test of random number 
generators [1,2]. Such a test has recently been performed by Wolff [14] for 
p = 223 and p = 389. In both cases no discrepancy between the simulation 
data and the exact analytic results was found. 

d. SU(2) lattice gauge theory. The generator with p = 223 is now being 
used in some high-precision calculations of the running coupling in the SU(2) 
lattice gauge theory [15]. So far all results obtained are compatible with earlier 
computations where shift register generators were employed. 

5.4 Recommended values of p 

From the theoretical discussion and the tests of the improved generator one 
concludes that the remaining statistical correlations are small when p is greater 
than about 200. The recommended default value is p = 223, and if one has 
any doubts that the simulation results might be biased by the random number 
generator, one may still set p = 389. A decorrelation of successive vectors of 
r random numbers down to the least significant digits is then guaranteed. 

To take still larger values of p appears to be pointless, since no empirical 
test or theoretical consideration indicates that a further improvement will be 
achieved. 
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Table 2. Average time needed to produce 1 new random number (p = 223) 



machine 



time [fj,s] 



SUN 10-41 
HP 9000/735 
CRAY YMP (1 CPU) 
APE-100 (1 node) 



5 
2 

0.7 

5 



5.5 Implementation and timing 

As discussed in sect. 2, the Marsaglia-Zaman algorithm can be implemented 
exactly using single precision floating point arithmetic. If random numbers 
between and 1 are desired, it is advantageous to work directly with the 
numbers x n /b instead of x n . No rounding is implied by this renormalization 
since b is a power of 2, i.e. the implementation remains exact. 

A portable FORTRAN code for the improved generator has been devel- 
oped by James [16] and is available through the CPC library. The name of the 
program is RANLUX. It comes with an initialization subroutine and further 
entry points to save and read the state of the generator. 

The generator has also been implemented on the APE-100 parallel com- 
puter [17]. The program may be obtained through anonymous ftp by dialing 
141 . 108 . 16 . 27 and copying the contents of the directory pub/random, or by 
writing to the author (luscher@ipsl02.desy.de). 

Since one uses only a fraction of the basic sequence of random numbers, 
the improved generator tends to be slow. For numerical simulations of lattice 
field theories, where large quantitites of random numbers are requested, it is 
hence important to take full advantage of any pipelining capabilities of the 
hardware. A difficulty here is that the Marsaglia-Zaman recursion (2.1), (2. 2) 
refers to the carry bit c n _i computed in the preceding step and so is not 
suitable for vectorization. 

The problem can be overcome by running several copies of the generator 
in parallel, with different initial values. The arithmetic operations are then 
pipelined horizontally, i.e. when looping over the copies. On the APE-100, for 
example, a good efficiency is achieved with 24 copies on each node. Some care 
should of course be paid to properly initialize the generators. In view of the 
astronomical period of the generator, the chances that any two of the copies 
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yield significantly correlated random numbers are however extremely slim. 

Some timing benchmarks for the improved generator with p = 223 are 
listed in table 2 [14,17]. The programs were written in FORTRAN and 
APESE, a high-level language for the APE-100. It is obvious that the numbers 
quoted depend on many technical details. They should hence be interpreted as 
a rough estimate of what can be achieved with a modest programming effort. 



6. Concluding remarks 

A well-known problem with random number generators is that their qual- 
ity is difficult to assess in any rigorous way. Some confidence in the reliability 
of any given generator can of course be gained by performing a large number 
of statistical tests. But doubts will always remain that the generator might 
fail in the next test. 

There exists an impressive list of classical dynamical systems which have 
been shown to be strongly chaotic. The states in these systems move ran- 
domly on time scales substantially greater than a certain characteristic time, 
related to the Liapunov exponent of the system. It should be emphasized that 
randomness can be given a precise mathematical meaning in this framework. 

The random number generator discussed in this paper may be considered 
a discrete approximation to such a chaotic dynamical system. A theoretical 
understanding of why the algorithm yields statistically independent random 
numbers is thus obtained. On longer time scales theoretical support for the 
good quality of the generator comes from the spectral test and the fact that the 
period can be shown to be extremely long. One might object that the generator 
is too slow for large scale applications. But other parts of the program are often 
much more costly so that the extra computer time needed for the generator 
is insignificant. One may also prefer to pay the price rather than taking any 
risk of producing corrupted data, especially when spending months of parallel 
computer time to a single project. 

I would like to thank Ulli Wolff for performing the Ising model tests and 
providing some of the timing benchmarks quoted in table 2. I am also indebted 
to Fred James for various useful informations and constant encouragement. 
Helpful discussions with Kari Kankaala, Rainer Sommer, Marcus Speh, Frank 
Steiner and Peter Weisz are gratefully acknowledged. 
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